source("script_02_create_objects_for_graphs.R")
library(foreign)
library(plyr)

#read in senator-level data
send <- read.dta("simulated_sen_level_data_from_r_new.dta", convert.underscore=T)
send <- send[order(send$nominee, send$simulation),]
head(send)
bailey.send <-subset(send, bailey.observation==1)

#drop 3 nominees who were not voted on
nomd <- subset(nomd, nominee!="miers" & nominee!="ginsburg_d" & nominee!="thornberry")	
bailey.nomd <- subset(bailey.nomd, nominee!="miers" & nominee!="ginsburg_d" & nominee!="thornberry")	
nom.level.stats <- subset(nom.level.stats, nominee!="miers" & nominee!="ginsburg_d" & nominee!="thornberry")	
bailey.nom.level.stats <- subset(bailey.nom.level.stats, nominee!="miers" & nominee!="ginsburg_d" & nominee!="thornberry")	
	
##################################################################################################################
##SENATE VOTING ERRORS
##################################################################################################################
#DW
##################################################################################################################
#a) first, individual senators
#########################################################
send$senator.nom <- paste(send$senator.id, send$nominee)
#COURT-OUTCOME MODEL
#first, get predictions for each sim
send$senators.court.outcome.pred.dw <- ifelse( abs(send$senator.sim.dw - send$j5new.sim.dw) <= abs(send$senator.sim.dw-send$j5old.sim.dw),1,0)
send$court.outcome.error.dw <- ifelse( send$senators.court.outcome.pred.dw!=send$vote,1,0)
#get mean error at senator level
senator.level.mean.court.outcome.error.dw <- tapply(send$court.outcome.error.dw, send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.court.outcome.error.dw  <- tapply(send$court.outcome.error.dw, send$nominee, mean )

#for a given senator, how often do predictions differ across simulations
tt <- tapply(send$senators.court.outcome.pred.dw, send$senator.nom,mean)
mean(tt <=.1 | tt >=.9)
mean(tt >=.4 & tt <=.6)
#POSITION-TAKING MODEL
#first, get predictions for each sim
send$senators.pos.taking.pred.dw <- ifelse( abs(send$senator.sim.dw - send$nom.proj.sim.dw) <= abs(send$senator.sim.dw-send$j5old.sim.dw),1,0)
send$pos.taking.error.dw <- ifelse( send$senators.pos.taking.pred.dw!=send$vote,1,0)
#get mean error at senator level
senator.level.mean.pos.taking.error.dw  <- tapply(send$pos.taking.error.dw, send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.pos.taking.error.dw  <- tapply(send$pos.taking.error.dw, send$nominee, mean )
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
send$left.cut <- ifelse( send$above.dw==1,	(send$j5old.sim.dw + send$j5new.sim.dw)/2, ifelse(send$above.dw==0, (send$nom.proj.sim.dw + send$j5old.sim.dw)/2, NA))
send$right.cut <- ifelse( send$above.dw==0,	(send$j5old.sim.dw + send$j5new.sim.dw)/2, ifelse(send$above.dw==1, (send$nom.proj.sim.dw + send$j5old.sim.dw)/2, NA))
send$senators.mixed.pred.yes.dw <- ifelse(  (send$above.dw ==1 & send$senator.sim.dw  >= send$right.cut) | (send$above.dw==0 & send$senator.sim.dw  <=  send$left.cut ),1,0)
send$senators.mixed.pred.no.dw <- ifelse(  (send$above.dw ==1 & send$senator.sim.dw  < send$left.cut) | (send$above.dw==0 & send$senator.sim.dw  > send$right.cut ),1,0)
send$senators.mixed.pred.ambiguous <- ifelse(  (send$senator.sim.dw > send$left.cut) & (send$senator.sim.dw < send$right.cut),1,0)
send$senators.mixed.pred.dw <- ifelse(send$senators.mixed.pred.yes.dw == 1, 1, ifelse(send$senators.mixed.pred.no.dw==1,0,NA))#to just use unambigous votes

################################################################################
#b) 	#Using the Senate median as the unit of analysis
################################################################################
#COURT-OUTCOME
nomd$sen.median.court.outcome.pred.dw <- ifelse( abs(nomd$sen.median.sim.dw-nomd$j5new.sim.dw) <= abs(nomd$sen.median.sim.dw-nomd$j5old.sim.dw),1,0)
nom.level.stats$sen.median.court.outcome.pred.dw.mean <- tapply(nomd$sen.median.court.outcome.pred.dw, nomd$nominee,mean)
nom.level.stats$sen.median.court.outcome.pred.dw.binary <- ifelse(nom.level.stats$sen.median.court.outcome.pred.dw.mean >=.5,1,0)
nom.level.stats$sen.median.court.outcome.pred.dw.binary.correct <- ifelse(nom.level.stats$sen.median.court.outcome.pred.dw.binary==nom.level.stats$confirmed,1,0)
#POSITION-TAKING
nomd$sen.median.pos.taking.pred.dw <- ifelse( abs(nomd$sen.median.sim.dw-nomd$nom.proj.sim.dw) <= abs(nomd$sen.median.sim.dw-nomd$j5old.sim.dw),1,0)
nom.level.stats$sen.median.pos.taking.pred.dw.mean <- tapply(nomd$sen.median.pos.taking.pred.dw, nomd$nominee,mean)
nom.level.stats$sen.median.pos.taking.pred.dw.binary <- ifelse(nom.level.stats$sen.median.pos.taking.pred.dw.mean >=.5,1,0)
nom.level.stats$sen.median.pos.taking.pred.dw.binary.correct <- ifelse(nom.level.stats$sen.median.pos.taking.pred.dw.binary==nom.level.stats$confirmed,1,0)
#MIXED
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
nomd$left.cut <- ifelse( nomd$above.dw==1,	(nomd$j5old.sim.dw + nomd$j5new.sim.dw)/2, ifelse(nomd$above.dw==0, (nomd$nom.proj.sim.dw + nomd$j5old.sim.dw)/2, NA))
nomd$right.cut <- ifelse( nomd$above.dw==0,	(nomd$j5old.sim.dw + nomd$j5new.sim.dw)/2, ifelse(nomd$above.dw==1, (nomd$nom.proj.sim.dw + nomd$j5old.sim.dw)/2, NA))
nomd$sen.median.mixed.pred.dw.yes <- ifelse(  (nomd$above.dw ==1 & nomd$sen.median.sim.dw  >= nomd$right.cut) | (nomd$above.dw==0 & nomd$sen.median.sim.dw  <=  nomd$left.cut ),1,0)
nomd$sen.median.mixed.pred.dw.no <- ifelse(  (nomd$above.dw ==1 & nomd$sen.median.sim.dw  < nomd$left.cut) | (nomd$above.dw==0 & nomd$sen.median.sim.dw  > nomd$right.cut ),1,0)
nomd$sen.median.mixed.pred.dw.ambiguous <- ifelse(  (nomd$sen.median.sim.dw > nomd$left.cut) & (nomd$sen.median.sim.dw < nomd$right.cut),1,0)
nomd$sen.median.mixed.pred.dw <- ifelse(nomd$sen.median.mixed.pred.dw.yes == 1, 1, ifelse(nomd$sen.median.mixed.pred.dw.no==1,0,NA))#to just use unambigous votes

####################
####################
#now get summaries to make  2 by 2 with uncertainty across simulations
temp.df.1.dw <- ddply(send, .(simulation), summarize
	, senators.court.outcome.pred.dw.no.actual.no = mean( senators.court.outcome.pred.dw ==0 & vote==0, na.rm=T )
	, senators.court.outcome.pred.dw.yes.actual.no = mean( senators.court.outcome.pred.dw ==1 & vote==0, na.rm=T )
	, senators.court.outcome.pred.dw.no.actual.yes = mean( senators.court.outcome.pred.dw ==0 & vote==1, na.rm=T )	
	, senators.court.outcome.pred.dw.yes.actual.yes = mean( senators.court.outcome.pred.dw ==1 & vote==1, na.rm=T )				
	, senators.pos.taking.pred.dw.no.actual.no = mean( senators.pos.taking.pred.dw ==0 & vote==0, na.rm=T )
	, senators.pos.taking.pred.dw.yes.actual.no = mean( senators.pos.taking.pred.dw ==1 & vote==0, na.rm=T )
	, senators.pos.taking.pred.dw.no.actual.yes = mean( senators.pos.taking.pred.dw ==0 & vote==1, na.rm=T )	
	, senators.pos.taking.pred.dw.yes.actual.yes = mean( senators.pos.taking.pred.dw ==1 & vote==1, na.rm=T )		
				)        
temp.df.2.dw <- ddply(nomd, .(simulation), summarize
	, sen.median.court.outcome.pred.dw.no.actual.no = mean( sen.median.court.outcome.pred.dw ==0 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.dw.yes.actual.no = mean( sen.median.court.outcome.pred.dw ==1 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.dw.no.actual.yes = mean( sen.median.court.outcome.pred.dw ==0 & confirmed==1, na.rm=T )	
	, sen.median.court.outcome.pred.dw.yes.actual.yes = mean( sen.median.court.outcome.pred.dw ==1 & confirmed==1, na.rm=T )				
	, sen.median.pos.taking.pred.dw.no.actual.no = mean( sen.median.pos.taking.pred.dw ==0 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.dw.yes.actual.no = mean( sen.median.pos.taking.pred.dw ==1 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.dw.no.actual.yes = mean( sen.median.pos.taking.pred.dw ==0 & confirmed==1, na.rm=T )	
	, sen.median.pos.taking.pred.dw.yes.actual.yes = mean( sen.median.pos.taking.pred.dw ==1 & confirmed==1, na.rm=T )		
			)      
#do mixed separately--drop ambigous votes
temp.df.3.dw <- ddply(subset(send, senators.mixed.pred.ambiguous==0), .(simulation), summarize
	, senators.mixed.pred.no.dw.actual.no = mean( senators.mixed.pred.dw ==0 & vote==0  , na.rm=T )
	, senators.mixed.pred.yes.dw.actual.no = mean( senators.mixed.pred.dw ==1 & vote==0 , na.rm=T )
	, senators.mixed.pred.no.dw.actual.yes = mean( senators.mixed.pred.dw ==0 & vote==1 , na.rm=T )	
	, senators.mixed.pred.yes.dw.actual.yes = mean( senators.mixed.pred.dw ==1 & vote==1, na.rm=T )		
			)  
temp.df.4 <- ddply(subset(nomd,sen.median.mixed.pred.dw.ambiguous==0), .(simulation), summarize
	, sen.median.mixed.pred.no.actual.no = mean( sen.median.mixed.pred.dw ==0 & confirmed==0   , na.rm=T )
	, sen.median.mixed.pred.yes.actual.no = mean( sen.median.mixed.pred.dw ==1 & confirmed==0 , na.rm=T )
	, sen.median.mixed.pred.no.actual.yes = mean( sen.median.mixed.pred.dw ==0 & confirmed==1 , na.rm=T )	
	, sen.median.mixed.pred.yes.actual.yes = mean( sen.median.mixed.pred.dw ==1 & confirmed==1 , na.rm=T )		
			)			
temp <- list(temp.df.1.dw, temp.df.2.dw,  temp.df.3.dw, temp.df.4)
sen.direct.tests.data.dw <- join_all(temp, by="simulation")

##################################################################################################################
#bly
##################################################################################################################
#a) first, individual senators
#########################################################
bailey.send$senator.nom <- paste(bailey.send$senator.id, bailey.send$nominee)
#COURT-OUTCOME
#first, get predictions for each sim
bailey.send$senators.court.outcome.pred.bly <- ifelse( abs(bailey.send$senator.sim.bly - bailey.send$j5new.sim.bly) <= abs(bailey.send$senator.sim.bly-bailey.send$j5old.sim.bly),1,0)
bailey.send$court.outcome.error.bly <- ifelse( bailey.send$senators.court.outcome.pred.bly!=bailey.send$vote,1,0)

#get mean error at senator level
senator.level.mean.court.outcome.error.bly <- tapply(bailey.send$court.outcome.error.bly, bailey.send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.court.outcome.error.bly  <- tapply(bailey.send$court.outcome.error.bly, bailey.send$nominee, mean )

#for a given senator, how often do predictions differ across simulations
tt <- tapply(bailey.send$senators.court.outcome.pred.bly, bailey.send$senator.nom,mean)
mean(tt <=.1 | tt >=.9, na.rm=T)
mean(tt >=.4 & tt <=.6, na.rm=T)
#POSITION-TAKING
#first, get predictions for each sim
bailey.send$senators.pos.taking.pred.bly <- ifelse( abs(bailey.send$senator.sim.bly - bailey.send$nom.proj.sim.bly) <= abs(bailey.send$senator.sim.bly-bailey.send$j5old.sim.bly),1,0)
bailey.send$pos.taking.error.bly <- ifelse( bailey.send$senators.pos.taking.pred.bly!=bailey.send$vote,1,0)
#get mean error at senator level
senator.level.mean.pos.taking.error.bly  <- tapply(bailey.send$pos.taking.error.bly, bailey.send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.pos.taking.error.bly  <- tapply(bailey.send$pos.taking.error.bly, bailey.send$nominee, mean )
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
bailey.send$left.cut <- ifelse( bailey.send$above.bly==1,	(bailey.send$j5old.sim.bly + bailey.send$j5new.sim.bly)/2, ifelse(bailey.send$above.bly==0, (bailey.send$nom.proj.sim.bly + bailey.send$j5old.sim.bly)/2, NA))
bailey.send$right.cut <- ifelse( bailey.send$above.bly==0,	(bailey.send$j5old.sim.bly + bailey.send$j5new.sim.bly)/2, ifelse(bailey.send$above.bly==1, (bailey.send$nom.proj.sim.bly + bailey.send$j5old.sim.bly)/2, NA))
bailey.send$senators.mixed.pred.yes.bly <- ifelse(  (bailey.send$above.bly ==1 & bailey.send$senator.sim.bly  >= bailey.send$right.cut) | (bailey.send$above.bly==0 & bailey.send$senator.sim.bly  <=  bailey.send$left.cut ),1,0)
bailey.send$senators.mixed.pred.no.bly <- ifelse(  (bailey.send$above.bly ==1 & bailey.send$senator.sim.bly  < bailey.send$left.cut) | (bailey.send$above.bly==0 & bailey.send$senator.sim.bly  > bailey.send$right.cut ),1,0)
bailey.send$senators.mixed.pred.ambiguous <- ifelse(  (bailey.send$senator.sim.bly > bailey.send$left.cut) & (bailey.send$senator.sim.bly < bailey.send$right.cut),1,0)
bailey.send$senators.mixed.pred.bly <- ifelse(bailey.send$senators.mixed.pred.yes.bly == 1, 1, ifelse(bailey.send$senators.mixed.pred.no.bly==1,0,NA))#to just use unambigous votes


################################################################################
#b) 	#Using the Senate median as the unit of analysis
################################################################################
#COURT-OUTCOME
bailey.nomd$sen.median.court.outcome.pred.bly <- ifelse( abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$j5new.sim.bly) <= abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$j5old.sim.bly),1,0)
bailey.nom.level.stats$sen.median.court.outcome.pred.bly.mean <- tapply(bailey.nomd$sen.median.court.outcome.pred.bly, bailey.nomd$nominee,mean)
bailey.nom.level.stats$sen.median.court.outcome.pred.bly.binary <- ifelse(bailey.nom.level.stats$sen.median.court.outcome.pred.bly.mean >=.5,1,0)
bailey.nom.level.stats$sen.median.court.outcome.pred.bly.binary.correct <- ifelse(bailey.nom.level.stats$sen.median.court.outcome.pred.bly.binary==bailey.nom.level.stats$confirmed,1,0)
#POSITION-TAKING
bailey.nomd$sen.median.pos.taking.pred.bly <- ifelse( abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$nom.proj.sim.bly) <= abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$j5old.sim.bly),1,0)
bailey.nom.level.stats$sen.median.pos.taking.pred.bly.mean <- tapply(bailey.nomd$sen.median.pos.taking.pred.bly, bailey.nomd$nominee,mean)
bailey.nom.level.stats$sen.median.pos.taking.pred.bly.binary <- ifelse(bailey.nom.level.stats$sen.median.pos.taking.pred.bly.mean >=.5,1,0)
bailey.nom.level.stats$sen.median.pos.taking.pred.bly.binary.correct <- ifelse(bailey.nom.level.stats$sen.median.pos.taking.pred.bly.binary==bailey.nom.level.stats$confirmed,1,0)
#MIXED
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
bailey.nomd$left.cut <- ifelse( bailey.nomd$above.bly==1,	(bailey.nomd$j5old.sim.bly + bailey.nomd$j5new.sim.bly)/2, ifelse(bailey.nomd$above.bly==0, (bailey.nomd$nom.proj.sim.bly + bailey.nomd$j5old.sim.bly)/2, NA))
bailey.nomd$right.cut <- ifelse( bailey.nomd$above.bly==0,	(bailey.nomd$j5old.sim.bly + bailey.nomd$j5new.sim.bly)/2, ifelse(bailey.nomd$above.bly==1, (bailey.nomd$nom.proj.sim.bly + bailey.nomd$j5old.sim.bly)/2, NA))
bailey.nomd$sen.median.mixed.pred.bly.yes <- ifelse(  (bailey.nomd$above.bly ==1 & bailey.nomd$sen.median.sim.bly  >= bailey.nomd$right.cut) | (bailey.nomd$above.bly==0 & bailey.nomd$sen.median.sim.bly  <=  bailey.nomd$left.cut ),1,0)
bailey.nomd$sen.median.mixed.pred.bly.no <- ifelse(  (bailey.nomd$above.bly ==1 & bailey.nomd$sen.median.sim.bly  < bailey.nomd$left.cut) | (bailey.nomd$above.bly==0 & bailey.nomd$sen.median.sim.bly  > bailey.nomd$right.cut ),1,0)
bailey.nomd$sen.median.mixed.pred.bly.ambiguous <- ifelse(  (bailey.nomd$sen.median.sim.bly > bailey.nomd$left.cut) & (bailey.nomd$sen.median.sim.bly < bailey.nomd$right.cut),1,0)
bailey.nomd$sen.median.mixed.pred.bly <- ifelse(bailey.nomd$sen.median.mixed.pred.bly.yes == 1, 1, ifelse(bailey.nomd$sen.median.mixed.pred.bly.no==1,0,NA))#to just use unambigous votes
####################
####################
#now get summaries to make  2 by 2 with uncertainty across simulations
temp.df.1.bly <- ddply(bailey.send, .(simulation), summarize
	, senators.court.outcome.pred.bly.no.actual.no = mean( senators.court.outcome.pred.bly ==0 & vote==0, na.rm=T )
	, senators.court.outcome.pred.bly.yes.actual.no = mean( senators.court.outcome.pred.bly ==1 & vote==0, na.rm=T )
	, senators.court.outcome.pred.bly.no.actual.yes = mean( senators.court.outcome.pred.bly ==0 & vote==1, na.rm=T )	
	, senators.court.outcome.pred.bly.yes.actual.yes = mean( senators.court.outcome.pred.bly ==1 & vote==1, na.rm=T )				
	, senators.pos.taking.pred.bly.no.actual.no = mean( senators.pos.taking.pred.bly ==0 & vote==0, na.rm=T )
	, senators.pos.taking.pred.bly.yes.actual.no = mean( senators.pos.taking.pred.bly ==1 & vote==0, na.rm=T )
	, senators.pos.taking.pred.bly.no.actual.yes = mean( senators.pos.taking.pred.bly ==0 & vote==1, na.rm=T )	
	, senators.pos.taking.pred.bly.yes.actual.yes = mean( senators.pos.taking.pred.bly ==1 & vote==1, na.rm=T )		
				)        
temp.df.2.bly <- ddply(bailey.nomd, .(simulation), summarize
	, sen.median.court.outcome.pred.bly.no.actual.no = mean( sen.median.court.outcome.pred.bly ==0 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.bly.yes.actual.no = mean( sen.median.court.outcome.pred.bly ==1 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.bly.no.actual.yes = mean( sen.median.court.outcome.pred.bly ==0 & confirmed==1, na.rm=T )	
	, sen.median.court.outcome.pred.bly.yes.actual.yes = mean( sen.median.court.outcome.pred.bly ==1 & confirmed==1, na.rm=T )				
	, sen.median.pos.taking.pred.bly.no.actual.no = mean( sen.median.pos.taking.pred.bly ==0 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.bly.yes.actual.no = mean( sen.median.pos.taking.pred.bly ==1 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.bly.no.actual.yes = mean( sen.median.pos.taking.pred.bly ==0 & confirmed==1, na.rm=T )	
	, sen.median.pos.taking.pred.bly.yes.actual.yes = mean( sen.median.pos.taking.pred.bly ==1 & confirmed==1, na.rm=T )		
			)      
#do mixed separately--drop ambigous votes
temp.df.3.bly <- ddply(subset(bailey.send, senators.mixed.pred.ambiguous==0), .(simulation), summarize
	, senators.mixed.pred.no.bly.actual.no = mean( senators.mixed.pred.bly ==0 & vote==0  , na.rm=T )
	, senators.mixed.pred.yes.bly.actual.no = mean( senators.mixed.pred.bly ==1 & vote==0 , na.rm=T )
	, senators.mixed.pred.no.bly.actual.yes = mean( senators.mixed.pred.bly ==0 & vote==1 , na.rm=T )	
	, senators.mixed.pred.yes.bly.actual.yes = mean( senators.mixed.pred.bly ==1 & vote==1, na.rm=T )		
			)  
temp.df.4 <- ddply(subset(bailey.nomd,sen.median.mixed.pred.bly.ambiguous==0), .(simulation), summarize
	, sen.median.mixed.pred.no.actual.no = mean( sen.median.mixed.pred.bly ==0 & confirmed==0   , na.rm=T )
	, sen.median.mixed.pred.yes.actual.no = mean( sen.median.mixed.pred.bly ==1 & confirmed==0 , na.rm=T )
	, sen.median.mixed.pred.no.actual.yes = mean( sen.median.mixed.pred.bly ==0 & confirmed==1 , na.rm=T )	
	, sen.median.mixed.pred.yes.actual.yes = mean( sen.median.mixed.pred.bly ==1 & confirmed==1 , na.rm=T )		
			)			
temp <- list(temp.df.1.bly, temp.df.2.bly,  temp.df.3.bly, temp.df.4)
sen.direct.tests.data.bly <- join_all(temp, by="simulation")

#########################
#THESE PRESENTS RESULTSS FOR TABLE 2
#########################
#now create output for making tables 
lo <- .025
hi <- .975
r <- 3
##roll calls
#DW
#COURT-OUTCOME
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(senators.court.outcome.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.court.outcome.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(senators.pos.taking.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.pos.taking.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#MIXED
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(senators.mixed.pred.no.dw.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.dw.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.dw.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.dw.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.dw.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.dw.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.mixed.pred.no.dw.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.dw.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.dw.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.dw.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.dw.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.dw.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#bailey
#COURT-OUTCOME
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(senators.court.outcome.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.court.outcome.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(senators.pos.taking.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.pos.taking.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#MIXED
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(senators.mixed.pred.no.bly.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.bly.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.bly.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.bly.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.bly.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.bly.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.mixed.pred.no.bly.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.bly.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.bly.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.bly.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.bly.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.bly.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))

#########################
#CONFIRMATION DECISIONS
##########################
#DW
#COURT-OUTCOME
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(sen.median.court.outcome.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.court.outcome.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.median.court.outcome.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.court.outcome.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(sen.median.pos.taking.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.pos.taking.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.median.pos.taking.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.pos.taking.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#BAILEY
########################
#COURT-OUTCOME
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(sen.median.court.outcome.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.court.outcome.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.median.court.outcome.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.court.outcome.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.median.court.outcome.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.court.outcome.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(sen.median.pos.taking.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.pos.taking.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.median.pos.taking.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.median.pos.taking.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.median.pos.taking.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.median.pos.taking.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))



###########################################################################
#ANALYSES IN APPENDIX
###########################################################################

###############################
#APPENDIX A.3.1 / TABLE A-2
#Section A.3.1 : The role of nominee quality and party in roll call votes
###############################

#Regress false negatives on individual quality and same party
#use senator-vote as unit, keeping only those where errors are consistent enough across simulations
#keep observations where the senator is predicted to vote no (under pure and quasi)
n.vars <- 4#number of predictors, plus constant
court.outcome.model.dw <- matrix(NA,nrow=n.sims, ncol=n.vars )
court.outcome.model.bly <- matrix(NA,nrow=n.sims, ncol=n.vars )
pos.taking.model.dw <- matrix(NA,nrow=n.sims, ncol=n.vars )
pos.taking.model.bly <- matrix(NA,nrow=n.sims, ncol=n.vars )
				
send <- within (send,{
	sen.j5old.dist <- abs(senator.sim.dw-j5old.sim.dw)
	sen.j5new.dist <- abs(senator.sim.dw-j5new.sim.dw)
	sen.nom.dist   <-  abs(senator.sim.dw-nom.proj.sim.dw)
})
bailey.send <- within (bailey.send,{
	sen.j5old.dist <- abs(senator.sim.dw-j5old.sim.dw)
	sen.j5new.dist <- abs(senator.sim.dw-j5new.sim.dw)
	sen.nom.dist   <-  abs(senator.sim.dw-nom.proj.sim.dw)
})

for (i in 1:n.sims){
	#COURT-OUTCOME /DW
	keep <- send$senators.court.outcome.pred.dw==0
	ok <- send$simulation==i 
	foo <- glm(vote ~ quality +  sameprty +  sen.nom.dist  
			, data=send,subset=ok & keep, family=binomial(link="logit"))
	court.outcome.model.dw[i,] <- coef(foo)
		#COURT-OUTCOME /bailey
	keep <- bailey.send$senators.court.outcome.pred.bly==0
	ok <- bailey.send$simulation==i 
	foo <- glm(vote ~ quality +  sameprty +  sen.nom.dist 
			, data=bailey.send,subset=ok & keep, family=binomial(link="logit"))
	court.outcome.model.bly[i,] <- coef(foo)
		#POSITION-TAKING /DW
	keep <- send$senators.pos.taking.pred.dw==0
	ok <- send$simulation==i 
	foo <- glm(vote ~ quality +  sameprty +   sen.nom.dist 
			, data=send,subset=ok & keep, family=binomial(link="logit"))
	pos.taking.model.dw[i,] <- coef(foo)
		#POSITION-TAKING /bailey
	keep <- bailey.send$senators.pos.taking.pred.bly==0
	ok <- bailey.send$simulation==i 
	foo <- glm(vote ~ quality +  sameprty +  +  sen.nom.dist  
			, data=bailey.send,subset=ok & keep, family=binomial(link="logit"))
	pos.taking.model.bly[i,] <- coef(foo)
}

#OUTPUT FOR REGRESSION TABLE
quantile(	court.outcome.model.dw[,1], c(.025,.5,.975))
quantile(	court.outcome.model.dw[,2], c(.025,.5,.975))
quantile(	court.outcome.model.dw[,3], c(.025,.5,.975))
quantile(	court.outcome.model.dw[,4], c(.025,.5,.975))

quantile(	court.outcome.model.bly[,1], c(.025,.5,.975))
quantile(	court.outcome.model.bly[,2], c(.025,.5,.975))
quantile(	court.outcome.model.bly[,3], c(.025,.5,.975))
quantile(	court.outcome.model.bly[,4], c(.025,.5,.975))

quantile(	pos.taking.model.dw[,1], c(.025,.5,.975))
quantile(	pos.taking.model.dw[,2], c(.025,.5,.975))
quantile(	pos.taking.model.dw[,3], c(.025,.5,.975))
quantile(	pos.taking.model.dw[,4], c(.025,.5,.975))

quantile(	pos.taking.model.bly[,1], c(.025,.5,.975))
quantile(	pos.taking.model.bly[,2], c(.025,.5,.975))
quantile(	pos.taking.model.bly[,3], c(.025,.5,.975))
quantile(	pos.taking.model.bly[,4], c(.025,.5,.975))

####################################################################
####################################################################
#A.4: EVALUATING MTM-THEORY OVER TIME
#graph probability of a mistake (by the Senate median) [DROPPING 3 nominees with no roll call votes]
	#if prob(confirm) > .5 and vote no, or pr(confirm)<.5 and vote yes, that's a mistake. =
#DROP 3 without roll call votes
temp <- subset(nom.level.stats, nominee!="thornberry" & nominee!="ginsburg_d"  & nominee!="miers" )
temp$court.outcome.prob.mistake.dw <- ifelse(temp$confirmed ==1, 1-temp$sen.median.court.outcome.pred.dw.mean, 
		temp$sen.median.court.outcome.pred.dw.mean )
temp$pos.taking.prob.mistake.dw <- ifelse(temp$confirmed ==1, 1-temp$sen.median.pos.taking.pred.dw.mean , 
		temp$sen.median.pos.taking.pred.dw.mean)

####now do errors at the senator level
#do proportion of votes that are predicted no that are actually yes
send$court.outcome.false.positive.dw <- ifelse( send$senators.court.outcome.pred.dw != send$vote,1,0)
temp$mean.court.outcome.false.negatives.dw <- tapply(send$court.outcome.error.dw, send$nominee, mean )
send$pos.taking.false.positive.dw <- ifelse( send$senators.pos.taking.pred.dw != send$vote,1,0)
temp$mean.pos.taking.false.negatives.dw <- tapply(send$pos.taking.error.dw, send$nominee, mean )

#combined graph of confirmation and roll call votes (FIGURE A-3)
pdf("Figure_A3.pdf", height=10, width=7)

banner.text <- 2.5
title.text <- 1.5
y.label.text <- 1.2
x.axis.text <- 1.0
y.axis.text <- 1.4
model.label.text <-2
bottom <- 4.5
left <- 4.5
top <- 0
right <- 0

axis.text <- function(axis=1,labels,at,srt=45,tcl=-0.4, cex.axis=1, col="black", font=1) {
 axis(axis,at=at,labels=FALSE,tcl=tcl, col=col)
 text(at,par("usr")[3]-.03,srt=srt,adj=1,labels=labels,xpd=TRUE, cex.axis=cex.axis, col=col, font=font)
}     

#SHOW REJECTIONS IN RED
nom.title.mistake <- temp$nominee.title.dw[order(temp$date.vote)]
nom.title.confirmed <- temp$confirmed[order(temp$date.vote)]

layout(cbind(c(1:6), c(1,7,8,4,9,10)), heights=c(.1,1,1,.01,1,1), widths=c(1,.05))
#CONFIRMATIONS
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(0,0, "Confirmation decisions", cex=banner.text, font=3, xpd=NA)
#box()

y <- temp$court.outcome.prob.mistake.dw[order(temp$date.vote)]
x <- 1:length(y)
offset <- .1
par( mar=c(bottom, left, top, right))
#COURT-OUTCOME
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red")

axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Pr(mistake) by Senate", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Pure Model", 3, line=0, font=2, cex=title.text)
lines(lowess( y~x, f=.75), col="dark green", lty=6)

#POSITION-TAKING
y <- temp$pos.taking.prob.mistake.dw[order(temp$date.vote)]
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red", font=1)

axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Pr(mistake) by Senate", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Quasi Model", 3, line=0, font=2, cex=title.text)
lines(lowess( y~x, f=.75), col="dark green", lty=6)

#ROLL CALLS
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(0,-5, "Individual roll call votes", cex=banner.text, font=3, xpd=NA)
#box()
y <-   temp$mean.court.outcome.false.negatives.dw[order(temp$date.vote)] #just false positives
offset <- .1
par( mar=c(bottom, left, top, right))
#COURT-OUTCOME
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
#ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red", font=1)

axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Prop. of false negatives", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Pure Model", 3, line=-1, font=2, cex=title.text)
lines(lowess(y~x, f=.5), col="dark green", lty=6)

y <-    temp$mean.pos.taking.false.negatives.dw[order(temp$date.vote)] #just false positives
#POSITION-TAKING
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
#ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red", font=1)
axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Prop. of false negatives", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Quasi Model", 3, line=-1, font=2, cex=title.text)
lines(lowess( y~x, f=.75), col="dark green", lty=6)

#label models
x.location <- -1
y.location <- -.3
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,y.location, "Court-outcome\nbased", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,y.location, "Position-taking\nsenators", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,y.location, "Court-outcome\nbased", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,y.location, "Position-taking\nsenators", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)

dev.off()

#Do for Bailey (not depicted in paper in the interests of space)
temp <- subset(bailey.nom.level.stats, nominee!="thornberry" & nominee!="ginsburg_d"  & nominee!="miers" )
temp$court.outcome.prob.mistake.bly <- ifelse(temp$confirmed ==1, 1-temp$sen.median.court.outcome.pred.bly.mean, 
		temp$sen.median.court.outcome.pred.bly.mean )
temp$pos.taking.prob.mistake.bly <- ifelse(temp$confirmed ==1, 1-temp$sen.median.pos.taking.pred.bly.mean , 
		temp$sen.median.pos.taking.pred.bly.mean)

		####now do errors at the senator level
#do proportion of votes that are predicted no that are actually yes
bailey.send$court.outcome.false.positive.bly <- ifelse(bailey.send$senators.court.outcome.pred.bly != bailey.send$vote,1,0)
temp$mean.court.outcome.false.negatives.bly <- tapply(bailey.send$court.outcome.error.bly, bailey.send$nominee, mean, na.rm=T )
bailey.send$pos.taking.false.positive.bly <- ifelse( bailey.send$senators.pos.taking.pred.bly != bailey.send$vote,1,0)
temp$mean.pos.taking.false.negatives.bly <- tapply(bailey.send$pos.taking.error.bly, bailey.send$nominee, mean, na.rm=T )

#This does the sae for Bailey (not depicted in paper in the interests of space)
pdf("Figure_A3_bailey_version.pdf", height=10, width=7)

banner.text <- 2.5
title.text <- 1.5
y.label.text <- 1.2
x.axis.text <- 1.0
y.axis.text <- 1.4
model.label.text <-2
bottom <- 4.5
left <- 4.5
top <- 0
right <- 0

axis.text <- function(axis=1,labels,at,srt=45,tcl=-0.4, cex.axis=1, col="black", font=1) {
 axis(axis,at=at,labels=FALSE,tcl=tcl, col=col)
 text(at,par("usr")[3]-.03,srt=srt,adj=1,labels=labels,xpd=TRUE, cex.axis=cex.axis, col=col, font=font)
}     

#SHOW REJECTIONS IN RED
nom.title.mistake <- temp$nominee.title.bly[order(temp$date.vote)]
nom.title.confirmed <- temp$confirmed[order(temp$date.vote)]

layout(cbind(c(1:6), c(1,7,8,4,9,10)), heights=c(.1,1,1,.01,1,1), widths=c(1,.01))
#CONFIRMATIONS
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(0,0, "Confirmation decisions", cex=banner.text, font=3, xpd=NA)
#box()

y <- temp$court.outcome.prob.mistake.bly[order(temp$date.vote)]
x <- 1:length(y)
offset <- .1
par( mar=c(bottom, left, top, right))
#COURT-OUTCOME
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red")

axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Pr(mistake) by Senate", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Pure Model", 3, line=0, font=2, cex=title.text)
lines(lowess( y~x, f=.75), col="dark green", lty=6)

#POSITION-TAKING
y <- temp$pos.taking.prob.mistake.bly[order(temp$date.vote)]
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red", font=1)

axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Pr(mistake) by Senate", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Quasi Model", 3, line=0, font=2, cex=title.text)
lines(lowess( y~x, f=.75), col="dark green", lty=6)

#ROLL CALLS
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(0,-5, "Individual roll call votes", cex=banner.text, font=3, xpd=NA)
#box()
y <-   temp$mean.court.outcome.false.negatives.bly[order(temp$date.vote)] #just false positives
offset <- .1
par( mar=c(bottom, left, top, right))
#COURT-OUTCOME
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
#ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red", font=1)

axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Prop. of false negatives", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Pure Model", 3, line=-1, font=2, cex=title.text)
lines(lowess(y~x, f=.5), col="dark green", lty=6)

y <-    temp$mean.pos.taking.false.negatives.bly[order(temp$date.vote)] #just false positives
#POSITION-TAKING
plot(x, y, type="n", axes=F, xlab="", ylab="", ylim=c(0,1))
segments(x-offset, 0,x-offset,y, lwd=2, col="black", lty=1)
#ok <- nom.title.confirmed==1
axis.text(1, at=x[ok], label=nom.title.mistake[ok], srt=45, cex.axis=x.axis.text, col="black")
axis.text(1, at=x[!ok], label=nom.title.mistake[!ok], srt=45, cex.axis=x.axis.text, col="red", font=1)
axis(2,mgp=c(2,.5,0), at=seq(0,1,.2), labels=c(0,".2",".4",".6",".8",1),las=1, cex.axis=y.axis.text)
mtext("Prop. of false negatives", 2,line=2.5, srt=90,cex=y.label.text)
#mtext("Quasi Model", 3, line=-1, font=2, cex=title.text)
lines(lowess( y~x, f=.75), col="dark green", lty=6)

#label models
x.location <- -4
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,0, "Median-driven", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,0, "Nominee-driven", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,0, "Median-driven", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)
par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
text(x.location,0, "Nominee-driven", cex=model.label.text, font=2, xpd=NA,srt=90, adj=0)

dev.off()

##################################################################################################################
#########REPEAT VOTING ERRORS AT SENATE LEVEL USING FILIBUSTER PIVOT
##################################################################################################################
#DW
################################################################################
#b) 	#Using the Senate median as the unit of analysis
################################################################################
nomd.pivots <- nomd
bailey.nomd.pivots <- bailey.nomd
nom.level.stats.pivots <- nom.level.stats
bailey.nom.level.stats.pivots <- bailey.nom.level.stats

#COURT-OUTCOME
nomd.pivots$sen.pivot.court.outcome.pred.dw <- ifelse( abs(nomd.pivots$sen.pivot.sim.dw-nomd.pivots$j5new.sim.dw) <= abs(nomd.pivots$sen.pivot.sim.dw-nomd.pivots$j5old.sim.dw),1,0)
nom.level.stats.pivots$sen.pivot.court.outcome.pred.dw.mean <- tapply(nomd.pivots$sen.pivot.court.outcome.pred.dw, nomd.pivots$nominee,mean)
nom.level.stats.pivots$sen.pivot.court.outcome.pred.dw.binary <- ifelse(nom.level.stats.pivots$sen.pivot.court.outcome.pred.dw.mean >=.5,1,0)
nom.level.stats.pivots$sen.pivot.court.outcome.pred.dw.binary.correct <- ifelse(nom.level.stats.pivots$sen.pivot.court.outcome.pred.dw.binary==nom.level.stats.pivots$confirmed,1,0)
#POSITION-TAKING
nomd.pivots$sen.pivot.pos.taking.pred.dw <- ifelse( abs(nomd.pivots$sen.pivot.sim.dw-nomd.pivots$nom.proj.sim.dw) <= abs(nomd.pivots$sen.pivot.sim.dw-nomd.pivots$j5old.sim.dw),1,0)
nom.level.stats.pivots$sen.pivot.pos.taking.pred.dw.mean <- tapply(nomd.pivots$sen.pivot.pos.taking.pred.dw, nomd.pivots$nominee,mean)
nom.level.stats.pivots$sen.pivot.pos.taking.pred.dw.binary <- ifelse(nom.level.stats.pivots$sen.pivot.pos.taking.pred.dw.mean >=.5,1,0)
nom.level.stats.pivots$sen.pivot.pos.taking.pred.dw.binary.correct <- ifelse(nom.level.stats.pivots$sen.pivot.pos.taking.pred.dw.binary==nom.level.stats.pivots$confirmed,1,0)
#MIXED
#define cutpoints (see Figure 2)
nomd.pivots$left.cut <- ifelse( nomd.pivots$above.dw==1,	(nomd.pivots$j5old.sim.dw + nomd.pivots$j5new.sim.dw)/2, ifelse(nomd.pivots$above.dw==0, (nomd.pivots$nom.proj.sim.dw + nomd.pivots$j5old.sim.dw)/2, NA))
nomd.pivots$right.cut <- ifelse( nomd.pivots$above.dw==0,	(nomd.pivots$j5old.sim.dw + nomd.pivots$j5new.sim.dw)/2, ifelse(nomd.pivots$above.dw==1, (nomd.pivots$nom.proj.sim.dw + nomd.pivots$j5old.sim.dw)/2, NA))
nomd.pivots$sen.pivot.mixed.pred.dw.yes <- ifelse(  (nomd.pivots$above.dw ==1 & nomd.pivots$sen.pivot.sim.dw  >= nomd.pivots$right.cut) | (nomd.pivots$above.dw==0 & nomd.pivots$sen.pivot.sim.dw  <=  nomd.pivots$left.cut ),1,0)
nomd.pivots$sen.pivot.mixed.pred.dw.no <- ifelse(  (nomd.pivots$above.dw ==1 & nomd.pivots$sen.pivot.sim.dw  < nomd.pivots$left.cut) | (nomd.pivots$above.dw==0 & nomd.pivots$sen.pivot.sim.dw  > nomd.pivots$right.cut ),1,0)
nomd.pivots$sen.pivot.mixed.pred.dw.ambiguous <- ifelse(  (nomd.pivots$sen.pivot.sim.dw > nomd.pivots$left.cut) & (nomd.pivots$sen.pivot.sim.dw < nomd.pivots$right.cut),1,0)
nomd.pivots$sen.pivot.mixed.pred.dw <- ifelse(nomd.pivots$sen.pivot.mixed.pred.dw.yes == 1, 1, ifelse(nomd.pivots$sen.pivot.mixed.pred.dw.no==1,0,NA))#to just use unambigous votes

####################
####################
#now get summaries to make  2 by 2 with uncertainty across simulations
temp.df.1.dw.using.pivots <- ddply(nomd.pivots, .(simulation), summarize
	, sen.pivot.court.outcome.pred.dw.no.actual.no = mean( sen.pivot.court.outcome.pred.dw ==0 & confirmed==0, na.rm=T )
	, sen.pivot.court.outcome.pred.dw.yes.actual.no = mean( sen.pivot.court.outcome.pred.dw ==1 & confirmed==0, na.rm=T )
	, sen.pivot.court.outcome.pred.dw.no.actual.yes = mean( sen.pivot.court.outcome.pred.dw ==0 & confirmed==1, na.rm=T )	
	, sen.pivot.court.outcome.pred.dw.yes.actual.yes = mean( sen.pivot.court.outcome.pred.dw ==1 & confirmed==1, na.rm=T )				
	, sen.pivot.pos.taking.pred.dw.no.actual.no = mean( sen.pivot.pos.taking.pred.dw ==0 & confirmed==0, na.rm=T )
	, sen.pivot.pos.taking.pred.dw.yes.actual.no = mean( sen.pivot.pos.taking.pred.dw ==1 & confirmed==0, na.rm=T )
	, sen.pivot.pos.taking.pred.dw.no.actual.yes = mean( sen.pivot.pos.taking.pred.dw ==0 & confirmed==1, na.rm=T )	
	, sen.pivot.pos.taking.pred.dw.yes.actual.yes = mean( sen.pivot.pos.taking.pred.dw ==1 & confirmed==1, na.rm=T )		
			)   
			   
#do mixed separately--drop ambigous votes
temp.df.2.dw.using.pivots <- ddply(subset(nomd.pivots,sen.pivot.mixed.pred.dw.ambiguous==0), .(simulation), summarize
	, sen.pivot.mixed.pred.no.actual.no = mean( sen.pivot.mixed.pred.dw ==0 & confirmed==0   , na.rm=T )
	, sen.pivot.mixed.pred.yes.actual.no = mean( sen.pivot.mixed.pred.dw ==1 & confirmed==0 , na.rm=T )
	, sen.pivot.mixed.pred.no.actual.yes = mean( sen.pivot.mixed.pred.dw ==0 & confirmed==1 , na.rm=T )	
	, sen.pivot.mixed.pred.yes.actual.yes = mean( sen.pivot.mixed.pred.dw ==1 & confirmed==1 , na.rm=T )		
			)			
temp <- list(temp.df.1.dw.using.pivots, temp.df.2.dw.using.pivots )
sen.direct.tests.data.dw.using.pivots <- join_all(temp, by="simulation")

##################################################################################################################
#BAILEY
##################################################################################################################
#COURT-OUTCOME
bailey.nomd.pivots$sen.pivot.court.outcome.pred.bly <- ifelse( abs(bailey.nomd.pivots$sen.pivot.sim.bly-bailey.nomd.pivots$j5new.sim.bly) <= abs(bailey.nomd.pivots$sen.pivot.sim.bly-bailey.nomd.pivots$j5old.sim.bly),1,0)
bailey.nom.level.stats.pivots$sen.pivot.court.outcome.pred.bly.mean <- tapply(bailey.nomd.pivots$sen.pivot.court.outcome.pred.bly, bailey.nomd.pivots$nominee,mean)
bailey.nom.level.stats.pivots$sen.pivot.court.outcome.pred.bly.binary <- ifelse(bailey.nom.level.stats.pivots$sen.pivot.court.outcome.pred.bly.mean >=.5,1,0)
bailey.nom.level.stats.pivots$sen.pivot.court.outcome.pred.bly.binary.correct <- ifelse(bailey.nom.level.stats.pivots$sen.pivot.court.outcome.pred.bly.binary==bailey.nom.level.stats.pivots$confirmed,1,0)
#POSITION-TAKING
bailey.nomd.pivots$sen.pivot.pos.taking.pred.bly <- ifelse( abs(bailey.nomd.pivots$sen.pivot.sim.bly-bailey.nomd.pivots$nom.proj.sim.bly) <= abs(bailey.nomd.pivots$sen.pivot.sim.bly-bailey.nomd.pivots$j5old.sim.bly),1,0)
bailey.nom.level.stats.pivots$sen.pivot.pos.taking.pred.bly.mean <- tapply(bailey.nomd.pivots$sen.pivot.pos.taking.pred.bly, bailey.nomd.pivots$nominee,mean)
bailey.nom.level.stats.pivots$sen.pivot.pos.taking.pred.bly.binary <- ifelse(bailey.nom.level.stats.pivots$sen.pivot.pos.taking.pred.bly.mean >=.5,1,0)
bailey.nom.level.stats.pivots$sen.pivot.pos.taking.pred.bly.binary.correct <- ifelse(bailey.nom.level.stats.pivots$sen.pivot.pos.taking.pred.bly.binary==bailey.nom.level.stats.pivots$confirmed,1,0)
#MIXED
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
bailey.nomd.pivots$left.cut <- ifelse( bailey.nomd.pivots$above.bly==1,	(bailey.nomd.pivots$j5old.sim.bly + bailey.nomd.pivots$j5new.sim.bly)/2, ifelse(bailey.nomd.pivots$above.bly==0, (bailey.nomd.pivots$nom.proj.sim.bly + bailey.nomd.pivots$j5old.sim.bly)/2, NA))
bailey.nomd.pivots$right.cut <- ifelse( bailey.nomd.pivots$above.bly==0,	(bailey.nomd.pivots$j5old.sim.bly + bailey.nomd.pivots$j5new.sim.bly)/2, ifelse(bailey.nomd.pivots$above.bly==1, (bailey.nomd.pivots$nom.proj.sim.bly + bailey.nomd.pivots$j5old.sim.bly)/2, NA))
bailey.nomd.pivots$sen.pivot.mixed.pred.bly.yes <- ifelse(  (bailey.nomd.pivots$above.bly ==1 & bailey.nomd.pivots$sen.pivot.sim.bly  >= bailey.nomd.pivots$right.cut) | (bailey.nomd.pivots$above.bly==0 & bailey.nomd.pivots$sen.pivot.sim.bly  <=  bailey.nomd.pivots$left.cut ),1,0)
bailey.nomd.pivots$sen.pivot.mixed.pred.bly.no <- ifelse(  (bailey.nomd.pivots$above.bly ==1 & bailey.nomd.pivots$sen.pivot.sim.bly  < bailey.nomd.pivots$left.cut) | (bailey.nomd.pivots$above.bly==0 & bailey.nomd.pivots$sen.pivot.sim.bly  > bailey.nomd.pivots$right.cut ),1,0)
bailey.nomd.pivots$sen.pivot.mixed.pred.bly.ambiguous <- ifelse(  (bailey.nomd.pivots$sen.pivot.sim.bly > bailey.nomd.pivots$left.cut) & (bailey.nomd.pivots$sen.pivot.sim.bly < bailey.nomd.pivots$right.cut),1,0)
bailey.nomd.pivots$sen.pivot.mixed.pred.bly <- ifelse(bailey.nomd.pivots$sen.pivot.mixed.pred.bly.yes == 1, 1, ifelse(bailey.nomd.pivots$sen.pivot.mixed.pred.bly.no==1,0,NA))#to just use unambigous votes

####################
####################
#now get summaries to make  2 by 2 with uncertainty across simulations
temp.df.1.bly.using.pivots <- ddply(bailey.nomd.pivots, .(simulation), summarize
	, sen.pivot.court.outcome.pred.bly.no.actual.no = mean( sen.pivot.court.outcome.pred.bly ==0 & confirmed==0, na.rm=T )
	, sen.pivot.court.outcome.pred.bly.yes.actual.no = mean( sen.pivot.court.outcome.pred.bly ==1 & confirmed==0, na.rm=T )
	, sen.pivot.court.outcome.pred.bly.no.actual.yes = mean( sen.pivot.court.outcome.pred.bly ==0 & confirmed==1, na.rm=T )	
	, sen.pivot.court.outcome.pred.bly.yes.actual.yes = mean( sen.pivot.court.outcome.pred.bly ==1 & confirmed==1, na.rm=T )				
	, sen.pivot.pos.taking.pred.bly.no.actual.no = mean( sen.pivot.pos.taking.pred.bly ==0 & confirmed==0, na.rm=T )
	, sen.pivot.pos.taking.pred.bly.yes.actual.no = mean( sen.pivot.pos.taking.pred.bly ==1 & confirmed==0, na.rm=T )
	, sen.pivot.pos.taking.pred.bly.no.actual.yes = mean( sen.pivot.pos.taking.pred.bly ==0 & confirmed==1, na.rm=T )	
	, sen.pivot.pos.taking.pred.bly.yes.actual.yes = mean( sen.pivot.pos.taking.pred.bly ==1 & confirmed==1, na.rm=T )		
			)      
#do mixed separately--drop ambigous votes
temp.df.2.bly.using.pivots <- ddply(subset(bailey.nomd.pivots,sen.pivot.mixed.pred.bly.ambiguous==0), .(simulation), summarize
	, sen.pivot.mixed.pred.no.actual.no = mean( sen.pivot.mixed.pred.bly ==0 & confirmed==0   , na.rm=T )
	, sen.pivot.mixed.pred.yes.actual.no = mean( sen.pivot.mixed.pred.bly ==1 & confirmed==0 , na.rm=T )
	, sen.pivot.mixed.pred.no.actual.yes = mean( sen.pivot.mixed.pred.bly ==0 & confirmed==1 , na.rm=T )	
	, sen.pivot.mixed.pred.yes.actual.yes = mean( sen.pivot.mixed.pred.bly ==1 & confirmed==1 , na.rm=T )		
			)			
temp <- list(temp.df.1.bly.using.pivots, temp.df.2.bly.using.pivots )
sen.direct.tests.data.bly.using.pivots <- join_all(temp, by="simulation")

#now create output for making tables in Open Office
lo <- .025
hi <- .975
r <- 3
#########################
#CONFIRMATION DECISIONS
##########################
#how often do predictions differ across median versus fili. pivot
with(nom.level.stats.pivots, table(sen.median.court.outcome.pred.dw.binary, sen.pivot.court.outcome.pred.dw.binary))
with(bailey.nom.level.stats.pivots, table(sen.median.court.outcome.pred.bly.binary, sen.pivot.court.outcome.pred.bly.binary))
with(nom.level.stats.pivots, nominee[sen.median.court.outcome.pred.dw.binary==1 &  sen.pivot.court.outcome.pred.dw.binary==0])
with(nom.level.stats.pivots, nominee[sen.median.court.outcome.pred.dw.binary==0 &  sen.pivot.court.outcome.pred.dw.binary==1])
with(bailey.nom.level.stats.pivots, nominee[sen.median.court.outcome.pred.bly.binary==1 &  sen.pivot.court.outcome.pred.bly.binary==0])
with(bailey.nom.level.stats.pivots, nominee[sen.median.court.outcome.pred.bly.binary==0 &  sen.pivot.court.outcome.pred.bly.binary==1])

with(nom.level.stats.pivots, table(sen.median.pos.taking.pred.dw.binary, sen.pivot.pos.taking.pred.dw.binary))
with(bailey.nom.level.stats.pivots, table(sen.median.pos.taking.pred.bly.binary, sen.pivot.pos.taking.pred.bly.binary))
with(nom.level.stats.pivots, nominee[sen.median.pos.taking.pred.dw.binary==1 &  sen.pivot.pos.taking.pred.dw.binary==0])
with(nom.level.stats.pivots, nominee[sen.median.pos.taking.pred.dw.binary==0 &  sen.pivot.pos.taking.pred.dw.binary==1])
with(bailey.nom.level.stats.pivots, nominee[sen.median.pos.taking.pred.bly.binary==1 &  sen.pivot.pos.taking.pred.bly.binary==0])
with(bailey.nom.level.stats.pivots, nominee[sen.median.pos.taking.pred.bly.binary==0 &  sen.pivot.pos.taking.pred.bly.binary==1])

########################
#TO CREATE TABLE A-3
##################
#DW
#COURT-OUTCOME
with(sen.direct.tests.data.dw.using.pivots,
matrix(data=c(paste("[", round(quantile(sen.pivot.court.outcome.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.court.outcome.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.pivot.court.outcome.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.court.outcome.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.dw.using.pivots,
matrix(data=c(paste("[", round(quantile(sen.pivot.pos.taking.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.pos.taking.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.pivot.pos.taking.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.pos.taking.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#BAILEY
#COURT-OUTCOME
with(sen.direct.tests.data.bly.using.pivots,
matrix(data=c(paste("[", round(quantile(sen.pivot.court.outcome.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.court.outcome.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.pivot.court.outcome.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.court.outcome.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.court.outcome.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.court.outcome.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.bly.using.pivots,
matrix(data=c(paste("[", round(quantile(sen.pivot.pos.taking.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.pos.taking.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(sen.pivot.pos.taking.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(sen.pivot.pos.taking.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(sen.pivot.pos.taking.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(sen.pivot.pos.taking.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))

#################################################
#ANALYSIS FOR APPENDIX SECTION A.6
		#REPLICATE SENATE VOTING DROPPING VOICE VOTES
##################################################################################################################
##SENATE VOTING ERRORS
##################################################################################################################
#DW
##################################################################################################################
#a) first, individual senators
#########################################################
#NB: I AM OVERWRITING DATA HERE, AND NOT RENAMING OBJECTS. SO YOU MUST RE-RUN FROM SCRATCH ABOVE TO GET RESULTS IN PAPER
send <- subset(send, !is.na(vote.adjusted))
bailey.send <- subset(bailey.send, !is.na(vote.adjusted))
nomd <- subset(nomd, rcvote==1)
bailey.nomd <- subset(bailey.nomd, rcvote==1)
nom.level.stats <- subset(nom.level.stats, rcvote==1)
bailey.nom.level.stats <- subset(bailey.nom.level.stats, rcvote==1)

send$senator.nom <- paste(send$senator.id, send$nominee)
#COURT-OUTCOME
#first, get predictions for each sim
send$senators.court.outcome.pred.dw <- ifelse( abs(send$senator.sim.dw - send$j5new.sim.dw) <= abs(send$senator.sim.dw-send$j5old.sim.dw),1,0)
send$court.outcome.error.dw <- ifelse( send$senators.court.outcome.pred.dw!=send$vote,1,0)
#get mean error at senator level
senator.level.mean.court.outcome.error.dw <- tapply(send$court.outcome.error.dw, send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.court.outcome.error.dw  <- tapply(send$court.outcome.error.dw, send$nominee, mean )

#POSITION-TAKING
#first, get predictions for each sim
send$senators.pos.taking.pred.dw <- ifelse( abs(send$senator.sim.dw - send$nom.proj.sim.dw) <= abs(send$senator.sim.dw-send$j5old.sim.dw),1,0)
send$pos.taking.error.dw <- ifelse( send$senators.pos.taking.pred.dw!=send$vote,1,0)
#get mean error at senator level
senator.level.mean.pos.taking.error.dw  <- tapply(send$pos.taking.error.dw, send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.pos.taking.error.dw  <- tapply(send$pos.taking.error.dw, send$nominee, mean )
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
send$left.cut <- ifelse( send$above.dw==1,	(send$j5old.sim.dw + send$j5new.sim.dw)/2, ifelse(send$above.dw==0, (send$nom.proj.sim.dw + send$j5old.sim.dw)/2, NA))
send$right.cut <- ifelse( send$above.dw==0,	(send$j5old.sim.dw + send$j5new.sim.dw)/2, ifelse(send$above.dw==1, (send$nom.proj.sim.dw + send$j5old.sim.dw)/2, NA))
send$senators.mixed.pred.yes.dw <- ifelse(  (send$above.dw ==1 & send$senator.sim.dw  >= send$right.cut) | (send$above.dw==0 & send$senator.sim.dw  <=  send$left.cut ),1,0)
send$senators.mixed.pred.no.dw <- ifelse(  (send$above.dw ==1 & send$senator.sim.dw  < send$left.cut) | (send$above.dw==0 & send$senator.sim.dw  > send$right.cut ),1,0)
send$senators.mixed.pred.ambiguous <- ifelse(  (send$senator.sim.dw > send$left.cut) & (send$senator.sim.dw < send$right.cut),1,0)
send$senators.mixed.pred.dw <- ifelse(send$senators.mixed.pred.yes.dw == 1, 1, ifelse(send$senators.mixed.pred.no.dw==1,0,NA))#to just use unambigous votes

#for a given senator, how often do predictions differ across simulations
tt <- tapply(send$senators.pos.taking.pred.dw, send$senator.nom,mean)
mean(tt <.4 | tt >.6)

################################################################################
#b) 	#Using the Senate median as the unit of analysis
################################################################################
#COURT-OUTCOME
nomd$sen.median.court.outcome.pred.dw <- ifelse( abs(nomd$sen.median.sim.dw-nomd$j5new.sim.dw) <= abs(nomd$sen.median.sim.dw-nomd$j5old.sim.dw),1,0)
nom.level.stats$sen.median.court.outcome.pred.dw.mean <- tapply(nomd$sen.median.court.outcome.pred.dw, nomd$nominee,mean)
nom.level.stats$sen.median.court.outcome.pred.dw.binary <- ifelse(nom.level.stats$sen.median.court.outcome.pred.dw.mean >=.5,1,0)
nom.level.stats$sen.median.court.outcome.pred.dw.binary.correct <- ifelse(nom.level.stats$sen.median.court.outcome.pred.dw.binary==nom.level.stats$confirmed,1,0)
#POSITION-TAKING
nomd$sen.median.pos.taking.pred.dw <- ifelse( abs(nomd$sen.median.sim.dw-nomd$nom.proj.sim.dw) <= abs(nomd$sen.median.sim.dw-nomd$j5old.sim.dw),1,0)
nom.level.stats$sen.median.pos.taking.pred.dw.mean <- tapply(nomd$sen.median.pos.taking.pred.dw, nomd$nominee,mean)
nom.level.stats$sen.median.pos.taking.pred.dw.binary <- ifelse(nom.level.stats$sen.median.pos.taking.pred.dw.mean >=.5,1,0)
nom.level.stats$sen.median.pos.taking.pred.dw.binary.correct <- ifelse(nom.level.stats$sen.median.pos.taking.pred.dw.binary==nom.level.stats$confirmed,1,0)
#MIXED
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
nomd$left.cut <- ifelse( nomd$above.dw==1,	(nomd$j5old.sim.dw + nomd$j5new.sim.dw)/2, ifelse(nomd$above.dw==0, (nomd$nom.proj.sim.dw + nomd$j5old.sim.dw)/2, NA))
nomd$right.cut <- ifelse( nomd$above.dw==0,	(nomd$j5old.sim.dw + nomd$j5new.sim.dw)/2, ifelse(nomd$above.dw==1, (nomd$nom.proj.sim.dw + nomd$j5old.sim.dw)/2, NA))
nomd$sen.median.mixed.pred.dw.yes <- ifelse(  (nomd$above.dw ==1 & nomd$sen.median.sim.dw  >= nomd$right.cut) | (nomd$above.dw==0 & nomd$sen.median.sim.dw  <=  nomd$left.cut ),1,0)
nomd$sen.median.mixed.pred.dw.no <- ifelse(  (nomd$above.dw ==1 & nomd$sen.median.sim.dw  < nomd$left.cut) | (nomd$above.dw==0 & nomd$sen.median.sim.dw  > nomd$right.cut ),1,0)
nomd$sen.median.mixed.pred.dw.ambiguous <- ifelse(  (nomd$sen.median.sim.dw > nomd$left.cut) & (nomd$sen.median.sim.dw < nomd$right.cut),1,0)
nomd$sen.median.mixed.pred.dw <- ifelse(nomd$sen.median.mixed.pred.dw.yes == 1, 1, ifelse(nomd$sen.median.mixed.pred.dw.no==1,0,NA))#to just use unambigous votes

tt <- tapply(nomd$sen.median.court.outcome.pred.dw , nomd$nominee,mean)
mean(tt <.1 | tt >.9)
####################
####################
#now get summaries to make  2 by 2 with uncertainty across simulations
temp.df.1.dw <- ddply(send, .(simulation), summarize
	, senators.court.outcome.pred.dw.no.actual.no = mean( senators.court.outcome.pred.dw ==0 & vote==0, na.rm=T )
	, senators.court.outcome.pred.dw.yes.actual.no = mean( senators.court.outcome.pred.dw ==1 & vote==0, na.rm=T )
	, senators.court.outcome.pred.dw.no.actual.yes = mean( senators.court.outcome.pred.dw ==0 & vote==1, na.rm=T )	
	, senators.court.outcome.pred.dw.yes.actual.yes = mean( senators.court.outcome.pred.dw ==1 & vote==1, na.rm=T )				
	, senators.pos.taking.pred.dw.no.actual.no = mean( senators.pos.taking.pred.dw ==0 & vote==0, na.rm=T )
	, senators.pos.taking.pred.dw.yes.actual.no = mean( senators.pos.taking.pred.dw ==1 & vote==0, na.rm=T )
	, senators.pos.taking.pred.dw.no.actual.yes = mean( senators.pos.taking.pred.dw ==0 & vote==1, na.rm=T )	
	, senators.pos.taking.pred.dw.yes.actual.yes = mean( senators.pos.taking.pred.dw ==1 & vote==1, na.rm=T )		
				)        
temp.df.2.dw <- ddply(nomd, .(simulation), summarize
	, sen.median.court.outcome.pred.dw.no.actual.no = mean( sen.median.court.outcome.pred.dw ==0 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.dw.yes.actual.no = mean( sen.median.court.outcome.pred.dw ==1 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.dw.no.actual.yes = mean( sen.median.court.outcome.pred.dw ==0 & confirmed==1, na.rm=T )	
	, sen.median.court.outcome.pred.dw.yes.actual.yes = mean( sen.median.court.outcome.pred.dw ==1 & confirmed==1, na.rm=T )				
	, sen.median.pos.taking.pred.dw.no.actual.no = mean( sen.median.pos.taking.pred.dw ==0 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.dw.yes.actual.no = mean( sen.median.pos.taking.pred.dw ==1 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.dw.no.actual.yes = mean( sen.median.pos.taking.pred.dw ==0 & confirmed==1, na.rm=T )	
	, sen.median.pos.taking.pred.dw.yes.actual.yes = mean( sen.median.pos.taking.pred.dw ==1 & confirmed==1, na.rm=T )		
			)      
#do mixed separately--drop ambigous votes
temp.df.3.dw <- ddply(subset(send, senators.mixed.pred.ambiguous==0), .(simulation), summarize
	, senators.mixed.pred.no.dw.actual.no = mean( senators.mixed.pred.dw ==0 & vote==0  , na.rm=T )
	, senators.mixed.pred.yes.dw.actual.no = mean( senators.mixed.pred.dw ==1 & vote==0 , na.rm=T )
	, senators.mixed.pred.no.dw.actual.yes = mean( senators.mixed.pred.dw ==0 & vote==1 , na.rm=T )	
	, senators.mixed.pred.yes.dw.actual.yes = mean( senators.mixed.pred.dw ==1 & vote==1, na.rm=T )		
			)  
temp.df.4 <- ddply(subset(nomd,sen.median.mixed.pred.dw.ambiguous==0), .(simulation), summarize
	, sen.median.mixed.pred.no.actual.no = mean( sen.median.mixed.pred.dw ==0 & confirmed==0   , na.rm=T )
	, sen.median.mixed.pred.yes.actual.no = mean( sen.median.mixed.pred.dw ==1 & confirmed==0 , na.rm=T )
	, sen.median.mixed.pred.no.actual.yes = mean( sen.median.mixed.pred.dw ==0 & confirmed==1 , na.rm=T )	
	, sen.median.mixed.pred.yes.actual.yes = mean( sen.median.mixed.pred.dw ==1 & confirmed==1 , na.rm=T )		
			)			
temp <- list(temp.df.1.dw, temp.df.2.dw,  temp.df.3.dw, temp.df.4)
sen.direct.tests.data.dw <- join_all(temp, by="simulation")

##################################################################################################################
#bly
##################################################################################################################
#a) first, individual senators
#########################################################
bailey.send$senator.nom <- paste(bailey.send$senator.id, bailey.send$nominee)
#COURT-OUTCOME
#first, get predictions for each sim
bailey.send$senators.court.outcome.pred.bly <- ifelse( abs(bailey.send$senator.sim.bly - bailey.send$j5new.sim.bly) <= abs(bailey.send$senator.sim.bly-bailey.send$j5old.sim.bly),1,0)
bailey.send$court.outcome.error.bly <- ifelse( bailey.send$senators.court.outcome.pred.bly!=bailey.send$vote,1,0)
#get mean error at senator level
senator.level.mean.court.outcome.error.bly <- tapply(bailey.send$court.outcome.error.bly, bailey.send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.court.outcome.error.bly  <- tapply(bailey.send$court.outcome.error.bly, bailey.send$nominee, mean )

#for a given senator, how often do predictions differ across simulations
tt <- tapply(bailey.send$senators.court.outcome.pred.bly, bailey.send$senator.nom,mean)
mean(tt <=.1 | tt >=.9, na.rm=T)
mean(tt >=.4 & tt <=.6, na.rm=T)
#POSITION-TAKING
#first, get predictions for each sim
bailey.send$senators.pos.taking.pred.bly <- ifelse( abs(bailey.send$senator.sim.bly - bailey.send$nom.proj.sim.bly) <= abs(bailey.send$senator.sim.bly-bailey.send$j5old.sim.bly),1,0)
bailey.send$pos.taking.error.bly <- ifelse( bailey.send$senators.pos.taking.pred.bly!=bailey.send$vote,1,0)
#get mean error at senator level
senator.level.mean.pos.taking.error.bly  <- tapply(bailey.send$pos.taking.error.bly, bailey.send$senator.nom, mean )
#now at nominee level (averaging over sims)
nominee.level.mean.pos.taking.error.bly  <- tapply(bailey.send$pos.taking.error.bly, bailey.send$nominee, mean )
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
bailey.send$left.cut <- ifelse( bailey.send$above.bly==1,	(bailey.send$j5old.sim.bly + bailey.send$j5new.sim.bly)/2, ifelse(bailey.send$above.bly==0, (bailey.send$nom.proj.sim.bly + bailey.send$j5old.sim.bly)/2, NA))
bailey.send$right.cut <- ifelse( bailey.send$above.bly==0,	(bailey.send$j5old.sim.bly + bailey.send$j5new.sim.bly)/2, ifelse(bailey.send$above.bly==1, (bailey.send$nom.proj.sim.bly + bailey.send$j5old.sim.bly)/2, NA))
bailey.send$senators.mixed.pred.yes.bly <- ifelse(  (bailey.send$above.bly ==1 & bailey.send$senator.sim.bly  >= bailey.send$right.cut) | (bailey.send$above.bly==0 & bailey.send$senator.sim.bly  <=  bailey.send$left.cut ),1,0)
bailey.send$senators.mixed.pred.no.bly <- ifelse(  (bailey.send$above.bly ==1 & bailey.send$senator.sim.bly  < bailey.send$left.cut) | (bailey.send$above.bly==0 & bailey.send$senator.sim.bly  > bailey.send$right.cut ),1,0)
bailey.send$senators.mixed.pred.ambiguous <- ifelse(  (bailey.send$senator.sim.bly > bailey.send$left.cut) & (bailey.send$senator.sim.bly < bailey.send$right.cut),1,0)
bailey.send$senators.mixed.pred.bly <- ifelse(bailey.send$senators.mixed.pred.yes.bly == 1, 1, ifelse(bailey.send$senators.mixed.pred.no.bly==1,0,NA))#to just use unambigous votes

#for a given senator, how often do predictions differ across simulations
tt <- tapply(bailey.send$senators.pos.taking.pred.bly, bailey.send$senator.nom,mean)
mean(tt <.4 | tt >.6, na.rm=T)

################################################################################
#b) 	#Using the Senate median as the unit of analysis
################################################################################
#COURT-OUTCOME
bailey.nomd$sen.median.court.outcome.pred.bly <- ifelse( abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$j5new.sim.bly) <= abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$j5old.sim.bly),1,0)
bailey.nom.level.stats$sen.median.court.outcome.pred.bly.mean <- tapply(bailey.nomd$sen.median.court.outcome.pred.bly, bailey.nomd$nominee,mean)
bailey.nom.level.stats$sen.median.court.outcome.pred.bly.binary <- ifelse(bailey.nom.level.stats$sen.median.court.outcome.pred.bly.mean >=.5,1,0)
bailey.nom.level.stats$sen.median.court.outcome.pred.bly.binary.correct <- ifelse(bailey.nom.level.stats$sen.median.court.outcome.pred.bly.binary==bailey.nom.level.stats$confirmed,1,0)
#POSITION-TAKING
bailey.nomd$sen.median.pos.taking.pred.bly <- ifelse( abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$nom.proj.sim.bly) <= abs(bailey.nomd$sen.median.sim.bly-bailey.nomd$j5old.sim.bly),1,0)
bailey.nom.level.stats$sen.median.pos.taking.pred.bly.mean <- tapply(bailey.nomd$sen.median.pos.taking.pred.bly, bailey.nomd$nominee,mean)
bailey.nom.level.stats$sen.median.pos.taking.pred.bly.binary <- ifelse(bailey.nom.level.stats$sen.median.pos.taking.pred.bly.mean >=.5,1,0)
bailey.nom.level.stats$sen.median.pos.taking.pred.bly.binary.correct <- ifelse(bailey.nom.level.stats$sen.median.pos.taking.pred.bly.binary==bailey.nom.level.stats$confirmed,1,0)
#MIXED
#define cutpoints (see Figure "Predicted votes in the Mixed MTM model")
bailey.nomd$left.cut <- ifelse( bailey.nomd$above.bly==1,	(bailey.nomd$j5old.sim.bly + bailey.nomd$j5new.sim.bly)/2, ifelse(bailey.nomd$above.bly==0, (bailey.nomd$nom.proj.sim.bly + bailey.nomd$j5old.sim.bly)/2, NA))
bailey.nomd$right.cut <- ifelse( bailey.nomd$above.bly==0,	(bailey.nomd$j5old.sim.bly + bailey.nomd$j5new.sim.bly)/2, ifelse(bailey.nomd$above.bly==1, (bailey.nomd$nom.proj.sim.bly + bailey.nomd$j5old.sim.bly)/2, NA))
bailey.nomd$sen.median.mixed.pred.bly.yes <- ifelse(  (bailey.nomd$above.bly ==1 & bailey.nomd$sen.median.sim.bly  >= bailey.nomd$right.cut) | (bailey.nomd$above.bly==0 & bailey.nomd$sen.median.sim.bly  <=  bailey.nomd$left.cut ),1,0)
bailey.nomd$sen.median.mixed.pred.bly.no <- ifelse(  (bailey.nomd$above.bly ==1 & bailey.nomd$sen.median.sim.bly  < bailey.nomd$left.cut) | (bailey.nomd$above.bly==0 & bailey.nomd$sen.median.sim.bly  > bailey.nomd$right.cut ),1,0)
bailey.nomd$sen.median.mixed.pred.bly.ambiguous <- ifelse(  (bailey.nomd$sen.median.sim.bly > bailey.nomd$left.cut) & (bailey.nomd$sen.median.sim.bly < bailey.nomd$right.cut),1,0)
bailey.nomd$sen.median.mixed.pred.bly <- ifelse(bailey.nomd$sen.median.mixed.pred.bly.yes == 1, 1, ifelse(bailey.nomd$sen.median.mixed.pred.bly.no==1,0,NA))#to just use unambigous votes

####################
####################
#now get summaries to make  2 by 2 with uncertainty across simulations
temp.df.1.bly <- ddply(bailey.send, .(simulation), summarize
	, senators.court.outcome.pred.bly.no.actual.no = mean( senators.court.outcome.pred.bly ==0 & vote==0, na.rm=T )
	, senators.court.outcome.pred.bly.yes.actual.no = mean( senators.court.outcome.pred.bly ==1 & vote==0, na.rm=T )
	, senators.court.outcome.pred.bly.no.actual.yes = mean( senators.court.outcome.pred.bly ==0 & vote==1, na.rm=T )	
	, senators.court.outcome.pred.bly.yes.actual.yes = mean( senators.court.outcome.pred.bly ==1 & vote==1, na.rm=T )				
	, senators.pos.taking.pred.bly.no.actual.no = mean( senators.pos.taking.pred.bly ==0 & vote==0, na.rm=T )
	, senators.pos.taking.pred.bly.yes.actual.no = mean( senators.pos.taking.pred.bly ==1 & vote==0, na.rm=T )
	, senators.pos.taking.pred.bly.no.actual.yes = mean( senators.pos.taking.pred.bly ==0 & vote==1, na.rm=T )	
	, senators.pos.taking.pred.bly.yes.actual.yes = mean( senators.pos.taking.pred.bly ==1 & vote==1, na.rm=T )		
				)        
temp.df.2.bly <- ddply(bailey.nomd, .(simulation), summarize
	, sen.median.court.outcome.pred.bly.no.actual.no = mean( sen.median.court.outcome.pred.bly ==0 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.bly.yes.actual.no = mean( sen.median.court.outcome.pred.bly ==1 & confirmed==0, na.rm=T )
	, sen.median.court.outcome.pred.bly.no.actual.yes = mean( sen.median.court.outcome.pred.bly ==0 & confirmed==1, na.rm=T )	
	, sen.median.court.outcome.pred.bly.yes.actual.yes = mean( sen.median.court.outcome.pred.bly ==1 & confirmed==1, na.rm=T )				
	, sen.median.pos.taking.pred.bly.no.actual.no = mean( sen.median.pos.taking.pred.bly ==0 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.bly.yes.actual.no = mean( sen.median.pos.taking.pred.bly ==1 & confirmed==0, na.rm=T )
	, sen.median.pos.taking.pred.bly.no.actual.yes = mean( sen.median.pos.taking.pred.bly ==0 & confirmed==1, na.rm=T )	
	, sen.median.pos.taking.pred.bly.yes.actual.yes = mean( sen.median.pos.taking.pred.bly ==1 & confirmed==1, na.rm=T )		
			)      
#do mixed separately--drop ambigous votes
temp.df.3.bly <- ddply(subset(bailey.send, senators.mixed.pred.ambiguous==0), .(simulation), summarize
	, senators.mixed.pred.no.bly.actual.no = mean( senators.mixed.pred.bly ==0 & vote==0  , na.rm=T )
	, senators.mixed.pred.yes.bly.actual.no = mean( senators.mixed.pred.bly ==1 & vote==0 , na.rm=T )
	, senators.mixed.pred.no.bly.actual.yes = mean( senators.mixed.pred.bly ==0 & vote==1 , na.rm=T )	
	, senators.mixed.pred.yes.bly.actual.yes = mean( senators.mixed.pred.bly ==1 & vote==1, na.rm=T )		
			)  
temp.df.4 <- ddply(subset(bailey.nomd,sen.median.mixed.pred.bly.ambiguous==0), .(simulation), summarize
	, sen.median.mixed.pred.no.actual.no = mean( sen.median.mixed.pred.bly ==0 & confirmed==0   , na.rm=T )
	, sen.median.mixed.pred.yes.actual.no = mean( sen.median.mixed.pred.bly ==1 & confirmed==0 , na.rm=T )
	, sen.median.mixed.pred.no.actual.yes = mean( sen.median.mixed.pred.bly ==0 & confirmed==1 , na.rm=T )	
	, sen.median.mixed.pred.yes.actual.yes = mean( sen.median.mixed.pred.bly ==1 & confirmed==1 , na.rm=T )		
			)			
temp <- list(temp.df.1.bly, temp.df.2.bly,  temp.df.3.bly, temp.df.4)
sen.direct.tests.data.bly <- join_all(temp, by="simulation")

#now create output for making tables in Open Office
lo <- .025
hi <- .975
r <- 3
##roll calls
#DW
#COURT-OUTCOME
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(senators.court.outcome.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.court.outcome.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(senators.pos.taking.pred.dw.no.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.no.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.dw.yes.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.pos.taking.pred.dw.no.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.dw.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.dw.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.dw.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#MIXED
with(sen.direct.tests.data.dw,
matrix(data=c(paste("[", round(quantile(senators.mixed.pred.no.dw.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.dw.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.dw.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.dw.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.dw.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.dw.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.mixed.pred.no.dw.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.dw.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.dw.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.dw.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.dw.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.dw.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#bailey
#COURT-OUTCOME
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(senators.court.outcome.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.court.outcome.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.court.outcome.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.court.outcome.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.court.outcome.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#POSITION-TAKING				
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(senators.pos.taking.pred.bly.no.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.no.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.no.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.bly.yes.actual.no, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.yes.actual.no, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.yes.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.pos.taking.pred.bly.no.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.no.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.no.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.pos.taking.pred.bly.yes.actual.yes, lo),r),"],", " ", round(quantile(senators.pos.taking.pred.bly.yes.actual.yes, .5),r)," ,", "[", round(quantile(senators.pos.taking.pred.bly.yes.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))
#MIXED
with(sen.direct.tests.data.bly,
matrix(data=c(paste("[", round(quantile(senators.mixed.pred.no.bly.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.bly.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.bly.actual.no, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.bly.actual.no, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.bly.actual.no, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.bly.actual.no, hi),r),"]", sep=""), 
paste("[", round(quantile(senators.mixed.pred.no.bly.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.no.bly.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.no.bly.actual.yes, hi),r),"]", sep="") ,
paste("[", round(quantile(senators.mixed.pred.yes.bly.actual.yes, lo),r),"],", " ", round(quantile(senators.mixed.pred.yes.bly.actual.yes, .5),r)," ,", "[", round(quantile(senators.mixed.pred.yes.bly.actual.yes, hi),r),"]", sep="") 
)	, nrow=2, ncol=2, byrow=T,dimnames = list(c("Vote no", "Vote yes"),     c("Predicted no", "Predicted yes"))))










